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Abstract 


The electrode of a PEM fuel cell is a porous medium generally made of carbon cloth or paper. Such a porous electrode has been widely modeled 
as a homogeneous porous medium with a constant permeability in the literature of PEM fuel cell. In fact, most of gas diffusion media are not 
homogeneous having non-isotropic permeability. In case of carbon cloth, the porous structure consists of carbon fiber tows, the bundles of carbon 
fiber, and void spaces among tows. The combinational effect of the void space and tow permeability results in the effective permeability of the 
porous electrode. In this work, the lattice Boltzmann method is applied to the simulation of the flow in the electrode of a PEM fuel cell. The electrode 
is modeled as void space and porous region which has certain permeability and the Stokes and Brinkman equations are solved in the flow field using 
the lattice Boltzmann model. The effective permeability of the porous medium is calculated and compared to an analytical calculation showing a 
good agreement. It has been shown that the permeability of porous medium is strongly dependant on the fiber tow orientation in three-dimensional 
simulations. The lattice Boltzmann method is an efficient and effective numerical scheme to analyze the flow in a complicated geometry such as 


the porous medium. 
© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


The demand for energy has been rising in all its forms while 
supplies of key fossil fuels have shown signs of decline dur- 
ing the last few decades. The consequent rising prices are likely 
to push us towards alternative energy sources. Fuel cells and 
hydrogen technology have been at the center of highly active 
research with viable early markets among the renewable energy 
alternatives. The polymer electrolyte membrane (PEM) fuel cell 
is now regarded as a promising alternative of internal combus- 
tion engine owing to its competitive power density with high 
efficiency as well as zero tailpipe emission [1]. 

A PEM cell consists of a membrane electrolyte assembly 
(MEA) sandwiched between two bipolar plates as shown in 
Fig. 1. The gas diffusion layer (GDL), catalyst layer and poly- 
mer electrolyte membrane are referred to as MEA where current 
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is produced. Fuel and oxidant are supplied to both sides of MEA 
through the flow channels on the bipolar plates producing elec- 
tron in the anode catalyst layer and water in the cathode catalyst 
layer. The GDL is a porous medium generally made of carbon 
cloth or paper and plays an essential role in fuel cells of permit- 
ting gas to be transported from the flow channel to the catalyst 
layer. The GDL must permit liquid water to be transported from 
the catalyst layer into the flow channels to remove the liquid 
water from the cell. When the liquid water is accumulated in the 
GDL the gas transport from the gas flow channel to the catalyst 
layer is hindered limiting the performance of a PEM fuel cell. 
Such water flooding is acute at the cathode side where the water 
is formed. The GDL also has to supply sufficient conductivity 
to shuttle the electrons between the catalyst layer and the bipo- 
lar plate. Normally, the conductivity of the GDL is inversely 
proportional to the thickness and porosity. Therefore, finding 
optimal values of the GDL parameters such as thickness, poros- 
ity, permeability and wetting characteristics have been important 
issues among PEM fuel cell researchers [2—6]. The pores can 
be hydrophobic with optimal value of polytetrafluoroethylene 
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Nomenclature 

a the distance from the center to the edge of the 
ellipse on the major axis 

b the distance from the center to the edge of the 
ellipse on the semi-major axis 

c ôx/ôt 

Cs sound speed 

ej velocity vector for the lattice Boltzmann model 

f,(x,t) distribution function 


fx, t) equilibrium distribution function 
Kef effective permeability 
Kiow permeability tensor 


Ktow tow permeability 


L length of computational domain 
m molecular weight 

P pressure 

R universal gas constant 

Sap Einstein summation convention 
T absolute temperature 

u velocity vector 

Ux velocity in x direction 

Uy velocity in y direction 

x position vector 


Greek letters 


B a parameter controlling the magnitude of the 
momentum sink 

ôt time step 

ôx the distance between abutting lattices 

E effective porosity 

Emin the ratio of void area to the area of cross section 


at the plane x=25 
Enom nominal porosity 


Etow porosity of fiber tows 

u viscosity of fluid 

he effective viscosity 

v kinematic viscosity 
density of fluid 

Ty relaxation time 

wi weight coefficient for the lattice Boltzmann model 

Subscripts 

in inlet 

max maximum value 

min minimum value 

out outlet 

x x direction 

y y direction 

Z z direction 


(PTFE) not to be congested with liquid water [2]. Water flow 
in GDL was experimentally investigated for various thicknesses 
and pore sizes in Benziger et al. [3]. Paganin et al. [2] observed 
a performance decrement at higher current densities when the 


GDL thickness is increased. The thickness of the GDL was opti- 
mized by a mathematical modeling in Inoue et al. [4] and by cell 
performance tests in Lee et al. [5] and Jordan et al. [6], with 
various GDL parameters such as porosity and thickness. A thin 
GDL with small porosity results in good electrical conductiv- 
ity, however, efficient mass transport requires large pores. In a 
small and long serpentine flow channel, reactants can directly 
cross to neighboring channels due to high pressure gradient 
and short path [7—9]. It has been found that such cross flow 
plays an important role for an effective removal of liquid water 
from the gas diffusion layer [10]. The amount of cross flow is 
strongly correlated with the thickness and permeability of the 
GDL [7,8]. 

Present study presents the application of lattice Boltzmann 
(LB) method to the micro-scale flow simulation of reactant in the 
GDL of a PEM fuel cell. The LB method has been accepted as a 
new computational tool for a variety of fluid transport phenom- 
ena [11,12]. It was applied to incompressible fluid flows [13,14], 
transport of passive scalars [15], miscible and immiscible flu- 
ids in complex geometries [16] and two-phase flow with phase 
change [17]. The kinetic nature of the LB method was also shown 
to be applicable to simulation of chemical reaction in micro- and 
meso-scopic flow [18] and electrokinetic transport phenomena 
[19]. We used the porous LB model developed by Spaid and Phe- 
lan [20] to simulate the reactant flow in porous electrode made 
of carbon cloth, which is a heterogeneous porous medium. The 
model is equivalent to solving a Stokes/Brinkman formulation 
to allow fluid transport through finite sized porous objects by 
treating them as an effective medium with a known permeabil- 
ity. Due to the dynamic nature of the LB calculations, the method 
has advantage of being able to capture the flow dynamics leading 
up to the steady solution of the Stokes and Brinkman equations. 
The electrode is modeled as void space and circular (cylinder in 
3-D) porous region which has a certain value of permeability. 
The effective permeability of the media from 2-D simulation 
is validated by comparing the LB result with a semi-analytical 
lubrication model described by Phelan and Wise [21]. The effect 
of fiber tow orientation on the effective permeability of the 
porous medium has been investigated using three-dimensional 
LB model within the wide ranges of tow permeability and 
porosity. 


2. Numerical section 
2.1. Governing equations 


Modeling micro-scale flow in fibrous porous media differs 
from the fluid mechanics problems in homogenous porous media 
in that there exist void spaces between fiber tows, as well as 
porous regions comprised of the tow itself. A conventional treat- 
ment of flow problems in which there exists a combination of 
void and porous regions involves application of the Stokes and 
Brinkman equations. Fluid flow in the open regions is governed 
by the Stokes equation given by 


uV-u=VP (1) 
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Fig. 1. Schematic of a polymer electrolyte fuel cell assembly. 


where u is viscosity, u is velocity vector and P is pressure. Inside 
the porous tows, the flow is modeled by the Brinkman equation, 


Laa 
Kiow 
where Kow is the permeability tensor for the fiber bundles, 
and He is the effective viscosity. The velocity and pressure are 
volume-averaged in the porous region. He is assumed to be the 
same as u since their relationship is not of primary importance 
for the present work. It is evident from Eq. (2) that the Brinkman 
equation reduces to Darcy’s law far from the void spaces when 
velocity gradients are small. At the interface between the open 
and porous regions, the Stokes/Brinkman formulation satisfies 
the proper boundary conditions such as the continuity of velocity 
and stress [20]. Eqs. (2) and (3) are coupled with the continuity 
equation given by 


MeV? (u) — (u) = V(P) (2) 


V-u=0 (3) 


where the velocity in Eq. (3) must be replaced by the volume- 
averaged velocity (u) within the fiber tows. 


(0, 0), 


i= 
c | cos 
ej = ( 2 


for i= 0 


i 
m, sin 


i —4— 1/2 
2c (cos Lie, sin 


2.2. Lattice Boltzmann method 


The LB method is based on a finite number of identical par- 
ticles that go through collision and propagation successively on 
prefixed paths in space. The following single-component lat- 
tice Boltzmann equation describes evolution of the distribution 


-1 
a): for i = 1,2,3,4 


function, f;(x,t), with the BGK collision term [22,23], 


f= FRO 


Ty 


Si + eiôr, t + ôr) — fiX, t) = (4) 
where e;’s are the discrete velocities, 5, the time step and T, is 
the relaxation time. fj“ represents the equilibrium distribution 
of f; given as 


2 
eq e-u (e-u) u-u 
edo, u) = w;0 |1 = 
fi (ow) = wip | + RF + ART ORT 
4/9 i=0 
Hae 1/9 i=1,2,3,4 (5) 
1/36 i=5,6,7,8 


where w;’s are the associated weight coefficients, R the universal 
gas constant and T is the absolute temperature. In Fig. 2(a), 
the velocity vectors, e;, for the two-dimensional 9-speed model 
(D2Q9) are shown to be 


i-—4-1/2 . 
—~ r], for i=5,6,7,8 

2 
where c = 6,/6;, and ôx is the distance between lattice points. The 
speed of the sound is c/./3 in the phase space and accordingly 
we have RT =c2/3. The macroscopic number density, o(x,/), and 
the velocity, u(x,t), of the fluid are obtained as 


8 
p= X mfi, (7) 
i=0 
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Fig. 2. (a) A microscopic picture of the woven carbon cloth gas diffusion media (white square represents the computational domain and white dotted circle represents 


the region of tow), (b) a schematic of 2-D numerical model for porous medium. 


8 
pu= > mfiei, (8) 
i=0 
where the molecular weight, m, is assumed to be unity in the 
present work. The kinematic viscosity is related to the relaxation 
time as 


(2t, — 1) 8&2 
v= A iy 
In this work, 6, and ô; are set to be unity for convenience, and 
thus, the fluid pressure is given in terms of the fluid density as 
p=p!3. 

The 3-D LBM is also based on LBGK model and solves Eq. 
(4) in a three-dimensional space. We use a three-dimensional 
15-speed model (D3Q15) for the present work [25,26] in which 
the velocity vectors, e;, i=0, 1,...,14, are the column vectors 
of the following matrix: 


(9) 


0 1 -1 0 0 0 0 1 -1 1 -1 1 
E=/]0 0 0 1 -1 O 0 1 -1 1 -1 -lI 
00 0 0 0 1 -1 1 -1 -1 1 1 


The density per node, p, and the macroscopic flow velocity, 
U= (Ux, Uy, Uz), are defined in terms of the particle distribution 
function by 


14 
=e (11) 
i=0 

14 
$ fier = pu (12) 
i=0 

The equilibrium can be chosen as [26] 

= o - zou -u, (13) 
f= se + Spei -u+ Lale - u)? — Zou -u, 
i=1,...,6: ClassI (14) 


' 64 24" 16" 48 i 


i=7,...,14: ClassII (15) 


where Classes I and II indicate the lattices at the face and at 
the corner of 3-D computational lattice structure, respectively, 
as shown in Fig. 2(b). Through Chapman-Enskog expansion, 
the Eqs. (4) and (5) lead to the Navier-Stokes equations near the 
incompressible limit [24]. They are given by continuity equation, 


0 

Raa Bg ee (16) 

ot 

and the momentum equation, 

iC Pua) + ðplpuaug) = —du(c? P) + Ag(2vSup); (17) 
-1 1 -t 

1 -1 1 (10) 
—1]1 -1 1 

where the Einstein summation convention is used. 


Sap =(dqug+Ogug)/2 is the strain-rate tensor. The sound 
speed is cs = ./3/8 for the 3-D LB model. 


2.3. Lattice Boltzmann model for a porous medium 


To recover the Brinkman equation using a lattice Boltzmann 
method, it is necessary to modify the standard equilibrium dis- 
tribution function to reduce the magnitude of the momentum, 
leaving the direction of momentum unchanged. This may be 
achieved by altering the velocity u(x,t) in the equilibrium distri- 
bution function by incorporating a forcing term given by 


TE(x, t) 
p(x, t) 


where U replaces u in Eq. (5) for the equilibrium distribution 
function, and the variable s(x,t) is either 0 or 1 depending on 


U = u(x, t) + s(x) 


(18) 
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weather a lattice site is located in a void or porous region, 
respectively. Eq. (18) has also been used to incorporate a body 
force such as fluid—solid interaction potential [27], fluid—fluid 
interaction potential [28] and electrostatic potential between 
electronically charged ions [19]. The form of the forcing term 
F(x,1) which is needed to recover the Brinkman equation is given 
by 


F(x, t) = —Bp(x, t)u(x, t), (19) 


where £ is a parameter controlling the magnitude of the momen- 
tum sink. From Eqs. (10) and (11), the equilibrium velocity in a 
Brinkman site (s(x,f) = 1.0) is defined as 


U = u(x, t)\(1 — Br). (20) 


In Spaid and Phelan [20], a proof is given that the above 
formulation removes a specific amount of momentum in the 
collision process, and hence recovers the Brinkman equation. 
Then, the LB equation in a Brinkman site is given for the steady 
state as 


1 
vV7u — pu = —VP. (21) 
p 


Eq. (21) reproduces the Brinkman equation if B=v/Kjow in 
which both the kinematic viscosity and the tow permeability 
are expressed in lattice units [21]. 


2.4. Boundary conditions 


The practical flow situation in a PEM fuel cell employs some 
forms of flow regulations at the inlet, outlet and other boundaries. 
Common types of regulations are constant-velocity, constant 
pressure and no-slip wall. For the D2Q9 model, we utilized the 
bounce-back rule to determine unknown part of particle distri- 
bution functions [25]. As an example, it is considered that the 
boundary is aligned with the x-direction with f4, f7, fg pointing 
into the wall as shown in Fig. 2(a). After streaming, fo, fi, f3, f4, 
fı, fg are known. Supposing that uy and uy are specified on the 
boundary Eq. (7) and (8) are used to determine fo, f5, fo and p 
as follows: 


fpot+fst+ fo=ep-(fot fit ft fat fi + fa), (22) 


fs — fo = pux — (fi — fa — fa + fs), (23) 
fot fs + fo = puy + (fa+ fa + fs), (24) 
1 
P=; [fo + fit f +2 fat f+ fs) (25) 
5 


To close the system of equations above, we assume the bounce 
back rule for the non-equilibrium part of the particle distribution 


normal to the boundary (in this case, fy — po = f4— 7, 
With f2 known fs and fe can be found, thus, 
2 
fa = fat Zpuy (26) 
1 1 1 
fs = fiz A= A) t Pe + g (27) 


S- fi) 3Ps + ipuy 28) 

The collision is also applied to the boundary nodes. No slip 
wall can be achieved applying velocities on the wall equals zero 
(ux = uy = 0) on the boundary. Similarly, when the density (pres- 
sure) is known on the boundary and the velocity is assumed to 
be normal to the boundary (uy =0), the unknowns (uy, f5, fg) can 
be determined from Eqs. (22)—(25) as: 


[Lfo+ fi + fs +2(fa+ fo + fal 


fo = fs 


ty =I (29) 
Pin 
2; 
fo = fat zay (30) 
Speen 31 
fs=fh fi fa) + guy (31) 
_ 1 1 32 
fo = fs z f) + gy (32) 


We followed Maier et al. [26] to impose boundary conditions 
for the 3-D LB simulations. To apply the constant pressure, the 
incoming populations through boundary are determined using 
the first-order extrapolation rule as: 


fix, 1) — fie + ei, t) = —[fj, t) — fie + ej, 0] (33) 


where e; is the outward vectors, which is the mirror image of e; at 
the boundary. This extrapolated populations are used to calculate 
a provisional density, p, from Eq. (11). The mass addition is used 
next to achieve the prescribed density, pin, at the boundary. The 
mass is added to the incoming Class II links are obtained as: 


1 
fia ) = fj D- FP (34) 


where p = p — pin. The set of incoming class II links satisfies the 
relation, e;-e, = 1, for indices i € class II. The no slip boundary 
condition can be enforced by using the bounce-back operation 
to set provisional values for the unknown populations, and then 
by redistributing these populations to achieve a prescribed tan- 
gential velocity at the wall. The bounce-back operation defines 
the density at the wall, and the mass redistribution defines the 
tangential momentum. The entire operation is executed after 
streaming and before collision. More details can be found in 
Maier et al. [26]. 


3. Results and discussion 


The enhancement of the reactant transportation in the GDL 
has been an important issue to improve the performance of a 
PEM fuel cell especially at high current density. The flow distri- 
bution in the GDL is also strongly affected by the cross flow that 
is also known to be effective on the liquid water removal from the 
GDL [7-10]. Therefore, the characteristic of reactant flow in the 
GDL would be one of the most important factors for the design 
and optimization of a PEM fuel cell. The GDL has been conven- 
tionally modeled as a homogeneous porous medium, as pointed 
out previously, however, the homogeneous model is not proper 
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(a) 
Fig. 3. Lattice structures of the 2-D and 3-D lattice Boltzmann model: (a) D2Q9 and (b) D3Q15. 


for the investigation of the flow characteristics in a heteroge- 
neous porous medium such as the GDL made of woven carbon 
cloth in a PEM fuel cell. The permeability of GDL depends on its 
fiber tow orientation and its non-isotropic permeability is critical 
on the flow and concentration distribution in a PEM fuel cell. 
For further detailed investigations, the LB method is applied to 
the flow simulation in the GDL of a PEM fuel cell to verify the 
effect of fiber tow orientation on the effective permeability of 


the porous medium. 


3.1. Two dimensional LB simulations 


In this section, the two dimensional LB model for the het- 
erogeneous porous medium is validated by comparing the LB 
simulation results with the lubrication model [21]. The GDL is 
modeled as consisting of void space and fiber tow area which has 
certain permeability. Eqs. (4)-(9) and (18)—-(20) are solved for a 
transverse flow over 2-D porous circle as shown in Fig. 3(b). For 
a general shape of an elliptical porous tow, the effective porosity 
of the medium is given by the lubrication model as: 

mx ab 


SS 2 
J Lil; 


(35) 


(1 — £tow) 


where a, b, Ly and Ly are defined in Fig. 4, and £tow is the porosity 
of the fiber tow. It is often useful to refer to a nominal porosity, 
which is based on the tow shape alone, defined as 


mx ab 


T , 36 
4 ae eo) 


Enom = 1 — 


Flow is driven in the y direction by a pressure difference 
between the boundaries y=—L, and y=Ly, while periodic 
boundary conditions are imposed at x=—L, and x= Ly. Then, 
regardless of the amount of the applied pressure, the effective 
permeability for the test geometry can be found as [21] 


2Ly/Lx 


Ke © - 37 
f Ty, WF dy (37) 


where 


409 


Carbon fiber 
Porous area 


O 
o0%50 
9986 


Void space 


(b) 


h? —1/2;2 ~1/2 
f= z tKiowletowg + (1+ Etow)h + Ktow A^ coth(Kigy 2] 


(38) 


Here, g and A are defined in Fig. 4. 

The porous tow is modeled as a circle (a= b = d/2) in a square 
computational domain (Ly = Ly = L/2) for the present 2-D simu- 
lation. Fifty-one and 44 grid points are assigned for the length 
(L) of computational domain and the diameter (d) of the porous 
circle, respectively. Simulations are conducted with the constant 
value of 6 =0.125 while the value of tow permeability, Ktow, is 
changed within the range of 0.1-0.0001. The same boundary 
conditions described for the lubrication model are applied; con- 
stant pressure at y= —Ly (inlet) and y = Ly (outlet) and periodic 
condition at x=—L, and x = Ly. 

Fig. 5(a) presents the velocity and pressure distributions in 
the porous medium when the tow permeability, Ktow, is 0.01. 
The velocity vectors are normalized by the maximum velocity 
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Fig. 4. Geometry for transverse flow in the y direction over porous tows with 
elliptical cross section, indicating parameters used in the lubrication model [21]. 
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Fig. 5. 2-D lattice Boltzmann simulation results: (a) pressure contour with velocity vectors and (b) pressure contour with streamlines; Ktow = 0.01. 


in the computational domain where the tow area is marked as a 
dashed circle. It is seen that the magnitude of velocity is sub- 
stantially smaller inside the tow area than outside the tow area 
due to the flow resistance caused by the low permeability of the 
tow area. The pressure contours are symmetrical profiles for the 
upstream and downstream divided by a straight contour that is 
normal to the mean flow direction (y-direction). The magnitude 
of velocity is the largest at the middle of periodic boundaries 
(x= —L,, x= Ly) where the pressure gradient is also the largest. 
Fig. 5(b) shows the streamlines with the pressure contours in 
the porous medium. Although the magnitude of velocity in the 
porous medium is very small the flow permeates into the porous 
tow area and the resultant streamlines are normal to the pres- 
sure contour. As the value of tow permeability increases the 
tow behaves like a solid obstacle, and thus, the flow permeation 
becomes negligible inside the tow area. 
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Fig. 6. Comparison of the effective permeability between the LB simulation (cir- 
cles) and the lubrication model (square) [21], as a function of the dimensionless 
tow permeability; the nominal porosity (€nom) of porous medium is 0.4. 


In Fig. 6, the resultant effective permeabilities of the porous 
medium are compared between the LB simulation and the lubri- 
cation model to prove the validity of the present LB model. The 
effective permeability is calculated using Darcy’s law for the LB 
simulation, 


_ Ktow i V(P) 
H 


(u) = (39) 
where (u) is the average velocity in the computational domain. 
The LB simulation shows a fair agreement with the lubrication 
model although the deviation becomes noticeable at small per- 
meability, 1074. It seems that the hexagonal LB model in Spaid 
and Phelan [20] shows a better agreement with the results of 
lubrication model. This may be attributed to the fact that the 
hexagonal grid (D2Q7) fits better for the circular porous tow 
than the rectangular grid (D2Q9) does in the present simula- 
tion. The LB simulation involves considerable numerical errors 
for the value of tow permeability smaller than 1074. This is 
attributed to that the relaxation time, t, is also decreased along 
with the value of tow permeability as shown in Table 1. 


3.2. Three-dimensional LB simulations; effective 
permeability versus tow orientations 


The architecture of the fiber preform is one of the most impor- 
tant factors which determine the anisotropy of textile porous 
medium. Much effort has been made to measure the perme- 


Table 1 

Values of the relaxation time T, for various tow permeability Ktow 
Ktow Ty 

0.1 1.325 
0.01 0.5825 
0.001 0.50825 
0.0001 0.500825 
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Fig. 7. Schematic of 3-D model for the two different tow directions: (a) Case A—flow across the tow direction and (b) Case B—flow along the tow direction. 


ability of various materials for a wide range of flow conditions 
[29-34] and such experimental measurements were correlated 
with the types of fiber preforms in Jackson and James [35]. 
Textile porous medium such as woven cloth can be modeled by 
capturing features of reciprocating patterns in the fiber preforms. 
The effect of different architecture and resulting fiber orientation 
on effective permeability of porous medium was studied analyt- 
ically by Yu and Lee [36] in which a simplified permeability 
model was developed from one-dimensional Stokes/Brinkman 
equation based on the proposed rectangular unit cell geometry. 

In this work, the porous LB model is applied to the 3-D sim- 
ulation of the flow in fibrous porous medium to investigate the 
effect of the tow orientation on the effective permeability. As 
shown in Fig. 7, the fiber tow is modeled as a circular cylin- 
der which has certain permeability in a cubical computational 
domain consisted of 51 x 51 x 51 lattice points. The tow has 
a constant length of 0.92Z (equivalent to 47 lattice points) and 
various diameters 0.5—0.9L (equivalent to 25—45 lattice points). 
The point of origin is set at the corner of the cubical computa- 
tional domain as shown in Fig. 7. No-slip wall is assigned as the 
boundary conditions for the top and bottom surfaces (z=0 and 
z=L) and periodic condition for the right and left walls (y=0 
and y= L). Constant pressure is applied at the plane x=0 and 
x=L to drive the mean flow motion in the x-direction. Two dif- 
ferent flow configurations are considered in this study; Fig. 7(a) 
represents the case when the fiber tow is normal to the mean 
flow direction (x-direction) while Fig. 7(b) when the fiber tow 
is parallel with the mean flow direction. Hereinafter, the flow 
configurations in Fig. 7(a) and (b) are referred to as Case A and 
Case B, respectively. 

Fig. 8 shows the streamlines and the contour plots for the 
magnitude of the x-component velocity (U,) at the two differ- 
ent planes (x = 25 and y=25) to present the flow structures. The 
magnitude of velocity is normalized using the maximum veloc- 
ity in the computational domain. On the plane x = 25 in Fig. 8(a), 
it is seen that the magnitude of velocity is the largest between 
the wall and tow boundary and negligible inside the tow region. 
The void area changes along the mean flow direction for Case A 


since the fiber tow is normal to the mean flow direction. Case A 
has symmetrical velocity profiles on the plane y = 25 indicating 
that the flow is accelerated and decelerated before and after the 
plane x = 25 with the noticeable flow permeation into the porous 
region. On the other hand, the magnitude of x-component veloc- 
ity is consistent for Case B due to the fixed void area along the 
mean flow direction. The flow permeation into the porous region 
is almost negligible in Case B. 

In Fig. 9, the effective permeability of the porous medium is 
obtained by the LB simulation and presented as a function of tow 
permeability for various tow diameters, 0.5L-0.9L. The effec- 
tive permeability increases as the tow permeability increases for 
both the Case A and B. It is larger for the smaller tow diameter 
with the larger rate of increment against the increasing tow per- 
meability. Although Case A and B present similar trends against 
the increasing tow diameter and permeability Case B has sub- 
stantially larger effective permeability. This is mainly caused by 
that Case B has a much smaller minimal void area than Case A as 
shown in Table 2. For instance, the nominal porosity on the plane 
x= 25 for Case A and B in Fig. 7 are 0.74 and 0.5, respectively. 
Furthermore, it seems that the effect of the drag by the no-slip 
wall is more significant in Case A. The combinational effects of 
the smaller void area and larger wall drag result in substantially 
lower effective permeability for Case A. Similar phenomenon 
was reported by an analytical analysis in Yu and Lee [36] and 
by experimental measurements in Jackson and James [35]. The 
effect of this phenomenon could be significant on the flow and 


Table 2 

Nominal porosity of 3-D porous medium (€nom, the ratio of the void space to 
total volume) and the ratio of the void area to the area of cross-section at the 
plane x=25 (Emin) 


diow Enom Emin (Case A) Emin (Case B) 
0.5L 0.82 0.54 0.80 
0.6L 0.74 0.45 0.72 
0.7L 0.65 0.35 0.62 
0.8L 0.54 0.26 0.50 
0.9L 0.41 0.17 0.36 
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On the plane x = 25 
Ux 
0.07 0.14 0.29 0.43 


(a) 


On the plane x = 25 


Ux 
0.07 0.14 0.29 0.43 


(b) 


On the plane y = 25 


0.57 0.71 086 1.0 


On the plane y = 25 


0.57 0.71 0.86 1.0 


Fig. 8. Contour plots for the x-component velocity (Uy) with streamlines at two cross-sectional planes; dtow =0.8L, Ktow =0.05: (a) Case A and (b) Case B; the 


magnitude of velocity is normalized using maximum velocity in the domain. 


concentration distributions in the GDL of a PEM fuel cell since 
the cross flow can be enhanced or reduced depending on the 
local tow orientation. The liquid water accumulation may be 
increased when Case A is close to the U-turn of flow channel 
where the velocity of cross flow is much lower than other areas 
[7,10]. 

Fig. 10 compares the effective permeability of the porous 
medium between Case A and B as a function of tow diameter 
when the tow permeability is 0.1 and 0.01. The effective per- 
meability decreases as the tow diameter increases for the same 
value of tow permeability. Both the effect of fiber orientation 
and tow permeability are more significant when the diameter of 


the tow is smaller. The effect of tow direction is dominant over 
the effective permeability when the tow diameter is small while 
the effect of permeability becomes more significant as the tow 
diameter decreases. 

Fig. 11 shows the effective permeability of the porous 
medium as a function of nominal porosity when the permeabil- 
ity, Ktow, is 0.1 and 0.01. According to Eq. (35), the nominal 
porosity of the porous medium decreases as the tow permeabil- 
ity increases and the resultant difference between Kjow =0.1 and 
0.01 is about 9%. The effective permeability is distinct for dif- 
ferent tow permeability (0.1 and 0.01) shown in Fig. 10 while 
the data with different tow permeability seems to collapse onto 
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Effective Permeability 


0.0 0.1 0.2 0.3 0.4 0.5 
Tow Permeability 


Effective Permeability 


0.0 0.1 0.2 0.3 0.4 0.5 
Tow Permeability 
(b) 
Fig. 9. Effective permeability of the porous medium as a function of tow per- 


meability for various tow diameters: (a) Case A and (b) Case B. 


a single line when plotted as a function of the nominal porosity, 
but the lines are different for the two fiber orientations of paral- 
lel or normal to the flow direction, as clearly shown in Fig. 11. 
This might suggest that the actual effective permeability for the 
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Fig. 10. Comparison of the effective permeability of the porous medium as a 
function of tow diameter. 
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Fig. 11. Comparison of the effective permeability of the porous medium with 
experimental measurements available in literature; (M) Carman [29], ($) Wig- 
gins et al. [30], (4s) Sullivan [31], (¥ ) Brown [32], (# ) Bergelin et al. [33], 
(®) Ingmanson et al. [34]; all these experimental measurements are also cited 
and available in Jackson and James [35]. 


fibrous materials will be bounded by these two lines if the fiber 
is inclined to the flow direction. As seen in Fig. 11, most of 
the experimental data points [29-33] are indeed falling between 
these two lines. 


3.3. Numerical stability 


The LB method provides a valid alternative to solving the 
Stokes and Brinkman equations for elliptical equations. Since 
the LB method is inherently dynamic with the Stokes and 
Brinkman equations it would be a reasonable choice if the partic- 
ular application requires the unsteady flow dynamics. However, 
if only the steady results are desired, standard methods such as 
finite difference or finite element calculations could be more 
efficient. It has been found that the present method is numer- 
ically unstable for parameter values v and Ktow such that 6B = 
v/Ktow > 2 as also indicated in Spaid and Phelan [20]. The 
origin of the numerical instability is the depletion of distribu- 
tions components, which point in the direction of the velocity, 
ultimately leading to negative values for these distributions and 
numerical instability. Therefore, the case of an infinite momen- 
tum sink equivalent to an object with zero permeability must 
be treated with other type of technique such as wall boundary 
conditions. 

The results of the present 2-D simulation covers the range of 
non-dimensional effective permeability from 1074 to 1, which 
is equivalent to 10~!* to 1078 m? in a physical domain if the 
characteristic length scale is assumed as 100 um. This fairly 
covers the range for the permeability of the actual GDL in a 
PEM fuel cell in many literatures [7—10]. However, the present 
3-D simulation does not cover this range with the current num- 
ber of lattice points. The modeling flow through heterogeneous 
material having very low permeability requires larger lattices 
due to the stability constraints on £. For the case of heteroge- 
neous porous media, additional stability issues may arise when 
considering flow through a near impermeable object at a low 
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nominal porosity. In this case, the magnitude of the velocity 
may be quite large in the void fluid regions, and could cause 
numerical instability. Careful attention must be paid to the mag- 
nitude of the driving force in order to keep the velocity within 
the proper bounds. 


4. Concluding remarks 


In this study, the porous LB model for the Stokes/Brinkman 
formulation has been applied to the flow simulation in a het- 
erogeneous porous medium such as the GDL made of woven 
carbon cloth in a PEM fuel cell. The GDL is modeled as a void 
space and fiber tow region, which has certain permeability. The 
LB model was validated by comparing the LB simulation results 
with an analytical calculation showing a good agreement for a 
wide range of tow permeability. The effective permeability of the 
porous medium is compared with analytical and experimental 
results in literature and the large variance of the effective perme- 
ability for the same porosity is explained in terms of the effect 
of tow permeability. The 3-D LB model is applied to investi- 
gate the effect of fiber tow orientation on the resultant effective 
permeability of a porous medium in a wide range of tow diam- 
eter and permeability. For the cases investigated, the effective 
permeability is substantially larger, as much as doubled, for the 
fiber orientation parallel with the mean flow direction, when 
compared to the case of the fibers normal to the mean flow 
direction. 

One of the primary advantages of LB method for simulat- 
ing fluid flows as compared to traditional numerical methods 
is their ability to robustly model interfaces between two or 
more fluids. The mechanism of liquid water accumulation in the 
GDL is of great interest for PEM fuel cell researchers since the 
phenomenon is one of the main mechanisms, which limit the per- 
formance of a PEM fuel cell. Conventional homogeneous porous 
model is not proper for the investigation of the liquid water accu- 
mulation in a heterogeneous porous medium due to the fact that 
the liquid water behavior inside the fiber tow is completely dif- 
ferent from that in the void space. The binary fluid model with the 
present Stokes/Brinkman formulation implemented would pro- 
vide a useful tool to investigate the mechanism of liquid water 
accumulation in the GDL of a PEM fuel cell. 
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